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ABSTRACT 

We extend the model of iKing et alJ (|2004) for variability in black hole accretion discs, by 
taking proper account of the thermal properties of the disc. Because the degree of variability 
in the iKing et alJ (|2004) model depends sensitively on the ratio of disc thickness to radius, 
H/R, it is important to follow the time-dependence of the local disc structure as the variability 
proceeds. In common with previous authors, we develop a one-zone model for the local disc 
structure. We agree that radial heat advection plays an important role in determining the inner 
disc struc ture, and also fin d limit-cycle behaviour. When the stochastic magnetic dynamo 
model of lKinget aU f2004) is added to these models, we find similar variability behaviour to 
before. 

We are now better placed to put physical constraints on model parameters. In particular, 
we find that in order to be consistent with the low degree of variability seen in the thermal 
disc component of black hole binaries, we need to limit the energy density of the poloidal 
field that can be produced by local dynamo cells in the disc to less than a few percent of the 
energy density of the dynamo field within the disc itself. 
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1 INTRODUCTION 

Accretion powered x-ray sources, both on the galactic (AGN) 
and stellar (X-ray binaries) scales, display significant aperi- 
odic variability, or flickerin g , on a broad range of timescales 
(see for example iMcHardM Il988l: iBelloni & HasingeJ 199C ; 
van der Klisl 1 1994 Markowitz et alJ 120031 : IVaughan et alJ ho03 : 
McClintock & Remillard 2004). The origin of this variability is not 
well understood. However. lUttlevI 12004) and llJttlev et alJ 120051 
have pointed out that the existence of a strong linear relationship 
between the amplitude of the X-ray variability and the amplitude of 
the x-ray flux (Uttlev & McHardv 2001) - the rms-flux relation - 
can be used as a diagnost ic to d istinguish between various models. 
For example Uttlev et al. 1 2005) emphasise that the standard simple 
shot-noise models, where the light curve is produced by the sum- 
mation of ran domly occurrin g sh ots, or flares, canno t explain such a 
relation. Both Uttlev (2004) and lUttlev et alJ l2005|) no te that such 
a relation can be produced naturally by the model of Lvubarskii 
J 19971) where the variability is produced by variations in the accre- 
tion rate occurring at different radii which propagate inwards, and 
so modulate the energy release in the central x-ray emitting region. 
The major problem with this model, as noted by Lvubarskii 1 1997), 
is the physical origin of the accretion rate varia tions. In o r der fo r 
the model to work, as also emphasised bv lChurazov et al] 1200 ll) . 
it is necessary for the timescale of the variations at each radius to 
be at least as long as the viscous timescale at that radius. The most 



obvious origin for fluctuations in an accretion disc are turbulent 
or hydro-magnetic, perhaps associated with the disc dynamo, and 
these take place on the local dynamical timescale (~ Q. \ where 
fl is the angular velocity), which is shorter than the local viscous 
timescale by a factor of ~ a(H/R) 2 , where a ( <1) is the vis cosity 
parameter, and H/R < 1 the disc opening angle iPrinale 1981). 

A solution to this problem, in the form of an explicit physical 
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model, has been recently proposed by King et al. 1 2004), based on 
ideas put forward bv lLivio et al J 12003). They suggest that the local 
dynamo processes in the disc can affect the accretion rate by driv- 
ing angular momentum loss in the form of an outflow (wind or jet). 
They model the dynamo as a small-scale stochastic phenomenon, 
operating on roughly the local dynamical timescale. They then pos- 
tulate that large-scale outflow can only occur when small-scale ran- 
dom processes in neighbouring disc annuli give rise by chance to a 
coherent large-scale field. This occurs on m uch longer timescales, 
by a factor of order 2 R/H iLivio et al 120 03 \ This also provides an 
explanation for why the dominant flickering frequencies are typi- 
cally much longer than the dynamical timescales at the centre of 
the disc where most of the energy is released. 

However, iKing et aT] 12004 took a very idealised model for 
the disc in that they assumed that the disc had a constant thick- 
ness ratio of H/R. This simplified the computations for two main 
reasons. First, local disc structure equations could be ignored, be- 
cause there was no need to compute the local value of disc thick- 
ness. Second, because the local disc dynamos were assumed to 
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be spatially independent on a radial scale of AR ~ H, there was 
no need to vary the number of independent dynamos as the disc 
thickness varied with time. In this paper we take the first steps to- 
wards remedying this deficiency. For the time being we concern 
ourselves with a standard (optically thick, geometrically thin) ac- 
cretion disc, and thus any conclusions we have will be relevant 
mainly to the the rmally dominated (TD) state of the X-ray tran- 
sients (e.g. McClintoc k & Remillardl2004l) . We leave the extension 
to the more interesting, and more variable (low/hard) state, con- 
taining both disc and corona, to futur e work. In Section [2] we ex- 
plain how we generalise the work of iKing et alJ 120041) by includ- 
ing equations to compute the local disc structure, and also explain 
how we vary the grid spacings to follow the radial size of the local 
disc dynamo cells in such a way as to minimise unwanted numer- 
ical mixing. In Sections [5] and |4| we investigate the models which 
result from our equations when there is a steady external accretion 
rate, and in Section[5]compare our results with previous work in the 
field. In Section|6|we add the variability according to the stochastic 
magnetic model of King et al. 1 2004). In Section^ we discuss the 
observational data on the variability of the thermal disc component, 
to which our models are relevant, and in Section|H|we show that for 
expected values of the model parameters our models are consistent 
with the data. We present brief conclusions in Section|5| 



2 METHODOLOGY 

In this Section we present the input physics for our basic accre- 
tion disc models. For a more detailed discussion, see for example 
lFranketalJl2002l) . 

The disc is treated as one-dimensional, in the sense that we re- 
solve the disc in radial direction only, dividing the disc into annuli, 
and just use the one-zone approximation in the vertical direction. 
Thus, all variables (e.g. mass density p, temperature T etc.) have 
one value at each radius R, and we make no further assumptions 
about the vertical disc structure. Thus, for example the surface den- 
sity X is given in terms of density p and disc semi-thickness H by 
the relation 



I = 2pH 



(1) 



As we discuss in Section [5] it is possible to make more detailed 
assumptions about the vertical disc structure. However, in view of 
the large number of uncertainties in the basic physical processes 
involved here, we regard such extra complications as unnecessary 
for our current purposes. 

2.1 Surface density 

The viscous evolutio n of the disc is governed jPringld 1 1 98 it 
iLivio & Pringlel 19921) by the equation of continuity 



<9£ 1 d 

— + (ZRu R ) = (2) 

dt RdR 

and by the equation of conservation of angular momentum 

— fetf 2 fi K ) + (Uiu R R 2 a K ) = . (3) 

dt V > RdR^ > 2nR dR 

Here S is the surface density, u R the radial velocity, M = 2nLRu R 
the accretion rate, Q K = -^GmJW the Keplerian rotation frequency 
and G the viscous torque. 

We should note two things here. First, although we are dis- 
cussing accretion discs around black holes in this paper, we use 



purely Newtonian gravity, around a point mass M, and truncate the 
disc at an appropriate radius R[ n = 6GM/c 2 . Here again we take 
the view that the wealth of uncertainties inherent in the physical 
processes involved, and the complicated nature of the subsequent 
behaviour, means that adding the complication of using the proper 
space-time geometry is not warranted at this stage. Second, by tak- 
ing the angular velocity of the disc material to be fix = -\/GM/R 3 
we have also neglected the effects of any radial pressure gradient in 
the disc. We shall show below that this approximation is justified 
for the disc models presented here. 

As is apparent from equation [3] angular momentum is trans- 
ported either by advection or viscous torques. For Keplerian rota- 
tion the viscous torque is given by 



G = -3nv t T.R 2 Cl K , 



(4) 



where v, is the kinematic viscosity. 

Equations an d can De combined lPringlell98lh to give 



3L 
~di 



3d_ 
RdR 



R 1/2 —(v t XR 1/2 ) 
dR y ' 



(5) 



We integrate this equation in time applying a zero-torque con- 
dition at the inner boundary (R m = 3R$, where R$ = 2GM/c 2 ) 
and thus set E s . =3Rs = 0. At the outer boundary we feed the disc 
with a constant external accretion rate M cxt . For advection we use 
a first-order donor cell upwind scheme. The integration conserves 
mass up to machine accuracy, accounting for mass supply at R out 
and loss at R[ n . 



2.2 Radial force balance and hydrostatic equilibrium 

We assume throughout that the disc is in hydrostatic equilibrium 
in the vertical direction (i.e. perpendicular to the disc plane). In 
the radial direction, the central gravitational pull is balanced by the 
centrifugal force 



, GM 

n 2 K R = —. 

R 2 



(6) 



The vertical structure is treated in the one-zone- 
approximation. Thus we take the vertical component of the 
gravitational force of the black hole to be balanced by the vertical 
pressure gradient, i.e. 



dP 

dz 



(7) 



where g- is the vertical gravity. In terms of our one-zone model we 
may write 



and 



dP 

dz 



-Q.IR-. 
K R 



(8) 



(9) 



Using this, together with EquationQwe obtain the relationship ex 
pressing vertical hydrodynamic equilibrium in the form 

„ GM . 



4pR 3 ' 



(10) 



2.3 Energy equation 



We consider the heat content q of one half of an elemental disc 
annulus of width AR at radius R, and height H. Then, from the 
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usual thermodynamic relations, and assuming hydrostatic equilib- 
rium (Equation llOt . we find (see Appendix lAl that 



dq 



de 



-L=e + u R —+APH + P 
at oR 



d(2nRHu R ) 
~dR ' 



(11) 



Here e is the internal energy of the semi-annulus, and the last two 
terms come from the 'PdV work done on the semi-annulus. A = 
2nRAR represents the surface area of a disc annulus of width AR at 
a given distance R from the black hole. The energy equation is then 



or). 



(12) 



The source terms on the RHS of the energy equation Jilt are 
given by the viscous dissipation of the disc per unit annulus area, 

. 9 GM 

Q = r*iF> (13) 

and radiative losses, Q~, which we take to be of the form 

QT = ^fr\ (14) 

where the optical depth r is given in terms of the opacity k r by 

T= l -Y.K K (p,T). (15) 

To ensure the consistency of the model we need disc to be optically 
thick in the vertical direction, i.e. t » 1. 

We have written the numerical scheme so that internal energy 
e is conserved to machine accuracy, apart from losses on the inner 
and outer boundary and local source terms (i.e. PdV). The bound- 
ary condition at the outer disc edge for the energy equation is taken 
to be divergence free, i.e. energy is being put in the last cell as it is 
lost to the next innermost cell. At the inner boundary energy is lost, 
assumed captured by the black hole. 



2.4.2 Viscosity 

We use the standard g-v iscosity prescription put forward by 
Shakura &SunvaevliT97^ and take the the (r</>)-component of the 
stress tensor to be ?,<j = —aP. Relating this t o the rate-of -strain ten- 
sor, t hrough the viscosity pv, (e.g. Bisnovatvi-Kogan & Lovelace! 
2001) so that t ri p = -jv t pQ. K we see that the kinematic viscosity is 
related to a by the equation 



V ' 3 a pCl K ' 



(20) 



This viscosity prescription (eq. 1201 . using Eqns. i 10t . Q 
and 1181 . can also be written in the form 



y t = -a^Jyc s H. 



(21) 



If we envisage the viscosity v, in physical terms as the product of a 
characteristic turbulent length and turbulent velocity, /, and V t , and 
our disc is assumed to be geometrically thin, we can associate the 
characteristic velocity and length scale with a fraction of the local 
scale height and of the local sound speed, respectively. In order 
that the turbulence remains subsonic we see that the factor 2/3a -\fy 
needs to be smaller than, or of order, unity. 

The viscous torque J4} using the viscosity prescription <20t 
can be written as 



G = -4naPHR 



(22) 



2.4.3 Opacity 



We take the Rosseland me an opacity, k r , tabulated by the OPAL 
opacity proiect fRogers & I gksiaj|l^9_^jK^lesiasj^ Rogerslll99r^ 
for the solar composition of Grevesse & Noels ( 1993). We use the 
X = 0.7 set of their opacity tables (X is the hydrogen mass fraction 
of the matter). Without further notice we fix the metallicity to be 
Z = 0.02 (solar metallicity). 



2.4 Material functions 

2.4.1 Equation of state 

The total pressure P is given by the sum of gas and radiation pres- 
sure (since t » 1) 



k B T 4cr , 
P = P— + —T\ 
pirip 3c 



(16) 



while the specific internal energy (J of the mixture of monatomic 
gas & radiation is 



2 pm v cp 
The sound speed c s can be calculated using 
P = ypc], 

where the adiabatic index y is given by 



d\o%P\ \6-\2pp-\pl 



dlogp/ s 12 -f ftp 



(17) 



(18) 



(19) 



and is a slowly varying function of /?/>, which ranges from | to | as 
Pp ranges from to 1, where /3 P is the ratio of gas pressure to total 
pressure. 



2.5 Numerical Grid and calculation procedure 

We use a one-dimensional grid in radius with N points ranging from 
R m = 3R$ to R oul . We solve the time-dependent equations J2j 
and lilt using the mass m and internal energy e of a disc semi- 
annulus as dependent variables. 

In order to allow simple application of the conce pts of vari- 
able d ynamo-driven angular momentum presented by iKing et alJ 
(2004), it is necessary for the widths (AR) of the disc annuli to be 
comparable to the local disc thickness H . Within the context of the 
one-zone approximation, this is a reasonable assumption, since in 
this case it makes little physical sense to try to attribute physical 
reality to attempts at resolving structures significantly smaller than 
the scale height H of the disc. 1 

While R iu and R out are kept constant, the number of grid-points 
is time-dependent and the number of grid points is adjusted accord- 
ing to the following prescription (see Fig.0. We remove/insert grid 
points if the local radial extent of the annulus is larger/smaller than 
2H and 0.5H, respectively. Grid points are inserted by conserv- 
ing disc mass, angular momentum and energy e. In typical runs, N 
ranges between a few dozens and a few hundred points. 



' lHameurv et alJ (l998) demonstrate that formal mathematical conver- 
gence of the numerical scheme may require a grid resolution much less 
than the disc thickness. We note here that mathematical convergence does 
not necessarily imply more accurate modelling of physical reality. 
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log R/(3 R s ) 

Figure 1. The refinement of the numerical grid. The Figure shows log H/R 
versus log R as it changes in a typical run. Every time the local scale height 
H exceeds(falls short of) twice(half) the local grid cell width, the resolution 
is reduced(increased). Due do the binarity requirement of the grid (see text) 
the same grid points are restored after repeated refinement/coarsening. 



The setup of the grid is implemented as follows. At the start of 
the calculation there are 2 points, at R m and R out . Then we calculate 
the scale height H(R ou[ ) and compare it with the interval R out - Rm- 
If R oui -R m > 2H(R out ), we put a new grid point at R new = ^R BUt R m 
and calculate the scale height H(R ncw ). This procedure is recur- 
sively applied on both intervals. The setup is complete if there is 
no need to refinement any more. 

From a numerical point of view it is beneficial in preventing 
numerical mixing to use a scheme which retains the structure of the 
grid cells fixed as far as possible. Thus, in order to maintain a binary 
structure of the grid (yielding the same coarse interval at the same 
place after several refinements), we assign each grid point with a 
number N, initially N(R m ) = and N(R out ) = 1. As a refinement 
occurs, the new grid point gets the mean number of the adjacent 
grid points. If this number is not an integer, all N are multiplied by 
2. The condition of grid coarsening is accompanied by the condi- 
tion for keeping the grid binary, i.e. we only remove a point if the 
N(R 0Ma ) - N(R mmme ) = N(R Temme ) - N(R m „ et ) and N(R mna ) has to 
be an integer multiple of N(R outa ) - N(R ml , el ). This binarity of the 
grid ensures that, for example, the magnetic field is never mixed 
numerically across the entire grid, i.e. subsequent refinement and 
coarsening restores the same grid cell again. 

A typical coarsening/refinement can be seen in Fig. □ The 
A(logi?) correlate very well with the H/R ratio since A(logi?) = 
log((fl + AR)/R) = log(l + AR/R) « AR/R = H/R as required. 
Every A log H/R = 0.3 dex, i.e. at log H/R » - 1 , - 1 .3 and - 1 .6 the 
resolution changes, as H de-/increases by a factor of 2. 



2.6 Advection 



It is o ften the case (see, for example, the discussion in Prin gle et alj 
1 19861) that the advective terms in the energy equation (i.e. those 
terms containing the radial velocity u R on the RHS of Equation ! 1 1 > 
are omitted. However, it was realised early that even i n the case of 
dwarf nova discs these ter ms can make a difference (Faulkner et alJ 
Il983l) . And more recentlv. lAbramowicz et alJll988t) drew attention 
to the fact that the radial advection of heat can play an important 
role in the local energy balance in the disc, especially at high ac- 
cretion rates close to the black hole, and can, in particular, have a 



stabilising effect on the disc in that region. To illustrate this, we note 
that the terms responsible for the radial advection of energy can be 
written in terms of the radial entropy gra dient as if they corre spond 
to a local additional heat loss term (e.g. [Muchotrzeb & Paczvnskil 
ll982tlBisnovatvi-Kogan & Lovelacel200ll) . in the form 



2ad 



(23) 



MT as 

2nRdR' 

The radial entropy gradient can be expanded in terms of pressure 
and radial density and temperature gradients. In doing so we obtain 



Gad 



M P 



(24) 



2nR 2 p' 

where the dimensionless measure of the strength of the effect,^, 
includes the radial derivatives. Usually, Xad(P) is a slowly varying 
function close to unity. However, although in the gas pressure dom- 
inated regime this value is indeed close to unity, in the radiation 
pressure dominated regime Xwi can be as large as 10, 

Using the expressions for the internal energy and pressure in 
Appendix lAl we can express^ in terms of the radial derivatives 
of p and T for a stationary disc 



,„ 21 \d\06 

12 Bp\ — 

2 P ) dlogR 



T d log p 

- + (4-3/?,=)—^. 



(25) 



3 STATIONARY SOLUTIONS AND STABILITY 

Before investigating the time dependent disc behaviour, and before 
adding in the complications of stochastically driven accretion, we 
first investigate steady disc solutions, and the conditions for stabil- 
ity. 



3.1 Stationary solutions 

In the stationary case, with accretion rate M = const., we have 
d/dt = and thus the conservation of angular momentum (Eqn.|3J 
becomes 



Mf = 37rv,E, 

conservation of energy i 1 1 i becomes 

Q + ~ Or - 2ad = o, 

and vertical hydrostatic equilibrium I lOt remains 

GM . 

P= , 

4pR 3 



(26) 



(27) 



(28) 



where the function f(R) = 1 - (R/R in )? ensures compliance with 
the inner boundary condition. For a given accretion rate, and at a 
particular radius, these three equations can be written as two equa- 
tions, depending on pressure P, density p and temperature T. The 
energy equation yields 



M 3 , , P 
and hydrostatic equilibrium yields 



l6no-T 4 aP 
3k rP Q k M/ 



0, 



p 3 - . A = o. 



(29) 



(30) 



167r 2 a 2 

Then for a given accretion rate, and at a particular radius, these two 
equations, together with the equation of state (Equation ll6t can be 
solved for P, p and T. 
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Figure 2. Local equilibrium solutions at R = 30Rs for a range of accre- 
tion rates. We consider solutions for different values of viscosity parameter 
a, the black hole mass M and the strength of radial advection, Xnd- Parts 
with negative M - S slope are thermally and viscously unstable. Note the 
stabilising effect of advection (when^d # 0). 



We solve these equations with a 2D nested-intervals method 
as described in Mayer & Duschl ( 2005). In general we consider the 
conditions at a particular radius R, for various values of the accre- 
tion rate M. For illustrative purposes we also consider the solutions 
for three different particular values of the strength of the adiabatic 
heat flux^ ad , and for various values of the viscosity parameter a. 

We present our results in terms of the Eddington luminosity 
and the corresponding Eddington accretion rate. In terms of the 
Eddington luminosity, L Edd = 4nGMm p c/crj we may define a cor- 
responding Eddington accretion rate by 



where r\ K is the efficiency. For a typical value of r\ K 
below) we find in numerical terms 



M Edd = 2.68 • 10- 



M 
10M E 



M Q yr 



(3D 
1/12 (see 

(32) 



3.2 S-Curves 

In Fig.|5|we plot the relationship between accretion rate M and a 
times disc surface density S for different values of radially constant 
a and for different values of x*& an d of the central mass M. 

These plots give us direct information about the stability of the 
steady-state discs. There are two types of instability which these 
discs can be subject to - v iscous (e.g. Liahtman & Eardlev ll974f) 
and thermal I Prinsle 1976). We show in Appendix|S|that the stabil- 
ity criterion in both cases is the same. The discs are unstable when 



dlogM 



< 0. 



(33) 



, fllogS , PT 

From Figure|2|we see that in the absence of heat advection the 
discs become unstable (i.e. the M(Z) curves have negative gradient) 
as the accretion rates approach th e Eddington limit. But, as pointed 
out by lAbramowicz et alj 1 19881) . advection of heat (i.e. non-zero 
Xni) provides stability. For all combinations of a and M considered, 
we see that within the physical range for^ ad there are always values 
of the accretion rate for which the discs are unstable, whereas for 




log E [g cm" ] 

Figure 3. Local equilibrium solution at R = 30Rs for M = 10 6 M Q , a = 1, 
Xad = 20 and different metallicities Z = 0, 0.001, 0.02(solar) and Z = 0.1. 
The hydrogen content is fixed at X = 0.7. 



higher black hole mass M, higher a and lower XaA me range of 
accretion rates we get instability beco mes larger. 

It is well known (see, for example. lLasotal20 1 ) that the insta- 
bility causes the disc to undergo limit-cycle behaviour. We investi- 
gate this further below. 



3.3 Dependency on black hole mass and metallicity 

For black hole masses typical of X-ray binaries (M » 10 M Q ), 
characteristic mid plane temperatures close to the black hole are of 
the order of 10 7 K. 

Since MQr K oc T 4 /t r (cf. eg. 127 \ in the absence of advection, 
we get in scaled units (R in Schwarzschild Radii, M in units of 
Eddington accretion rate) and roughly constant optical depth that 
T oc M -1 ' 4 . For a black hole mass of 10 6 Mo we then expect the 
mid plane temperature of the disc to be of the order of a few times 
10 5 A". 

For 10 6 Mq the disc temperature reaches T ss 10 5 4 K for 
reasonable accretion rates. Then the di sc shows an a dditional in- 
stability caused by the "Z-Bump" (cf. ISeaton et alJll994h in the 
opacity for a metallicity larger than Z = 0.02. This so call ed Z- 
Bump is also responsi ble for pulsations in j3 Cep stars iSimonll982t 
Kiriakidis et al. 1992). We show a S-curve presenting this instabil- 
ity in Fig. [3] It is evident that the disc for this instability locally can 
jump from accretion rates of 10~ 4 . . . 10"" 3 to 10 _1 MEdd- 

Generally, the disc becomes more unstable with increasing 
mass of the black hole, i.e. the lower turning point in the "S-curve" 
moves to lower accretion rates. 



3.4 H/R in the presence of advection 

For the approximations used in this paper to be valid we require 
that the disc thickness be small compared to the radius. We note 
here that the addition of a local heat sink in the guise of advection 
implies that the disc is typically thinner than it would be if all the 
energy generated locally were radiated at the same radius. In the 
stationary case, when advection of heat strongly dominates radia- 
tive losses (so that <2 ad » (2 _ ), Equation 1291 becomes 



l{a^Rff-- x ^ 
4 p 



0. 



(34) 
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Then using hydrostatic equilibrium (Equation llO> we find that 



H 
R 



2L 



(35) 



Thus, for example, taking Xm = 10 as a representative value, 
and with f = \ (corresponding to R « lORs) we find H/R « 0.2. 
Indeed, we see that as long as^- ad > 4/3, H/R is smaller than unity. 

3.5 Keplerian rotation 

For similar reasons, if advective heat flow is important, the disc 
angular velocity stays close to Keplerian. If we include the radial 
pressure gradient in the radial momentum balance equation, we find 
that the disc angular velocity SI is given by 



Q 2 



GM PR 



where 



-6> 



dlogP 
Slog if' 



This may be written as 



p |^ +(4 -3^) fll0gr 



(36) 



(37) 



(38) 



dlogR r dlogR 

The dominant term in the radiation pressure dominated do- 
main (p P = 0) is, of course, the temperature gradient, and a repre- 
sentative value is £/> » 2. 

In 1361 we treat the extra term in the brackets as small devi- 
ation. Thus we can replace GM/R by (£1 K R) 2 and with the hydro- 
static equilibrium i lOi we get 



CI 1 



(H\ 2 



(39) 



When advective heat transport is strongly dominant we can use 
Equation <35t to deduce that 



a 1 



4*ad 



(40) 



(41) 



For small deviation (Q. - Q.k <k Q.k) this gives approximately 
Q l 3/fr 

Thus, using representative values, the deviation from Keplerian ro- 
tation when advective heat transport is strong never exceeds around 
4 per cent. The assumption of Keplerian rotation is therefore a rea- 
sonable one. 



4 TIME-DEPENDENCE - RESULTS 

We have seen above that even in the absence of stochastic mag- 
netic phenomena our accretion disc models are expected to display 
time-dependent behaviour for some ranges of the parameters. We 
have carried out a number of simulations for different values of a, 
different external (or mean) accretion rates and different black hole 
masses. 

In this section, for reasons of computational time, we use a 
logarithmically equidistant grid with N = 500 grid points without 
refinement and coarsening. We solve the equations with an implicit 
solver. Comparison with the standard explicit solver and the refine- 
ment and coarsening as described in Sect. l2.5l does not show signifi- 
cant deviations. However the computing time decreased drastically 
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Figure 4. M versus a showing the systems with limit cycles (to the lower 
right of each line) and without (upper left) for a 10 M G black hole. The 
line is only a rough indication since the grid in (M, a-) is very coarse. The 
points, for which actually calculations have been done, are indicated by the 
crosses. 



and thus enabled us to carry out a parameter study concerning the 
global stability of these discs with respect to the instabilities dis- 
cussed above. 



4.1 General behaviour 

Our time-dependent models depend on the mass of the black hole, 
the external accretion rate and the value of the viscosity parameter 
a. 

As initial conditions we set up a stationary disc, assuming 
Q + = Q~ and neglecting advection. Depending on the accretion 
rate and the viscosity parameter, the disc either continues to stay 
in this stationary state, or, if advection of heat is important, the in- 
ner disc re-adjusts itself to allow for the advection. As we see from 
Figure [5] if the viscosity is high enough, we expect limit cycles to 
appear with the inner disc oscillating between a hot state with high 
accretion rate and a cool state with lower accretion rate (see Sec- 
tionl3~2t. 

If the inner disc is such that the mean accretion rate produces 
instability, the inner disc spends its time trying to jump between the 
two stable branches of the S-curve. While it is on the upper branch, 
with advection important, the disc is steadily depleted. The surface 
density decreases while the accretion rate is higher than average. 
As the inner disc is depleted, some matter is transported outward in 
a heating front. The inner disc is radiation pressure dominated and 
the accretion rate there is radially constant. Matter then piles up 
ahead of this front and subsequently the inner disc begins to starve 
and the accretion rate decreases. When the front reaches the stable, 
gas pressure dominated region of the disc, the front slows down and 
the inner disc cools, becoming more gas pressure dominated, as the 
supply of matter is interrupted. The dissipated energy is adverted 
outwards and the accretion rate in the inner disc decreases, vary- 
ing radially as M cc R 2 . Then inner disc reheats, the accretion rate 
increases inwards, the pileup is accreted and eventually the cycle 
restarts with the inner disc again becoming depleted. 
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4.2 Influence of a and the black hole mass on the stability 

The parameters for the grid of models with constant a discs we ran 
are shown in Figure|4]for a 10 Mo black hole. Also shown in Fig- 
ure are those parameter values for which limit-cycle behaviour 
is found. The higher the viscosity parameter, the more unstable the 
disc becomes since then advection can only stabilise for higher ac- 
cretion rates. 

For higher black hole masses, there are several complica- 
tions. First, the radiation pr essure dominated zone and hence the 
Liehtman & Ear dlevl 1 19741) unstable region is larger. For 10 7 Mq 
this zone is so large (more than 1000 Schwarzschild radii), that 
there is likely to be an interaction between this instability and the 
usual dwarf nova instability, caused by the opacit y drop at hy- 
drogen ionisation. Thi s has yet to be explored (cf. Clarke| [l988t 
IClarke & Shieldsll989h . 

Second, as discussed previously, there is a "Z-Bump" insta- 
bility present which makes the results extremely dependent on the 
chosen metal content for the opacity. 

Finally the computational effort becomes increasingly larger. 
The disc in the lLightman & Eardle\J 1 19741) unstable region oscil- 
lates between a thick disc and a more and more ultra-thin disc 
(H/R x 10~ 4 . . . 10~ 3 ) which puts strong limits on the time-step. 

These restrictions led us to concentrate on the 10 M Q black 
hole case and to leave the higher mass case for further investi- 
gations. We plot a sample lightcurve and a power density spec- 
trum for a 10 6 M n black hole in F ig. Q The two different in- 
stabilities (Lieht man & Eardlevl 1974L and the "Z-Bump" instabil- 
ity ) operate on different tim escales and different amplitudes. The 
Liehtman & Eardlevl 1 19741) instability leads to the huge outbursts 
in luminosity, while the small fluctuations are due to the "Z-Bump" 
instability. For the example chosen this instability affects only the 
region around 100 Schwarzschild radii. The instability then only 
appears as small fluctuations in the total disc luminosity. Note that 
for computational feasibility we set the number of grid points arti- 
ficially low. In order to fulfil our resolution criterion (AR ss H, cf. 
Sect. 12. 5> . we would have needed N ~ 3000 points to properly re- 
solve the inner disc in the low-M state. Thus the lightcurve is only 
indicative of the complications that can occur. 

4.3 Effect of advection on the observed accretion rate 

By definition the radiative luminosity of the disc is 

Lrad = 2n f Q RdR. (42) 

JR m 

For an infinitely extended, steady standard disc (R out — > oo, no ra- 
dial advection of heat), the energy dissipated in the disc by vis- 
cous processes (using our inner 'black hole' radius R ir[ = 3Rs = 
6GM/c 2 ) is 

1 GMM 1 . , 
Lam = ~~ 7i = -7^ Mr. (43) 

Thus the efficiency ^ K of our Newtonian, standard black hole ac- 
cretion disc in Keplerian rotation extending from R iu = 3R S to in- 
finity in converting rest mass energy into radiation is t] K = 

In the absence of advective heat flow, for an assumed value of 
the efficiency r] K a measurement of the observed radiative flux, i.e. 
irad, yields an estimate of the accretion rate by setting L a & = L d iss- 
However, when advection of heat is important, we expect devia- 
tions from the standard L mA ec M relation. In particular, since some 
of the dissipated heat is now advected directly into the black hole, 




W [M E ddI 

Figure 5. The time-averaged 'observed' accretion rate, Mobs, deduced from 
the radiative luminosity L lal j, is plotted versus the external (actual) accretion 
rate M for a 10 M Q black hole for different values of a. For high actual 
accretion rates, some of the accretion energy is advected into the black hole. 



we expect the observed radiated flux to give an underestimate of the 
actual accretion rate. From the numerical simulations, we calculate 
the radiative luminosity L ra d(r) at a given time, and deduce from an 
'observed' accretion rate given by 

Mobs(f) = ^ , (44) 
or the time-average of this 

i r T 

(M obs ) = == LUt) dt . (45) 

In Fig. |3 we plot < M obs > versus the actual accretion rate M, both 
in units of the Eddington accretion rate. If the disc model at the 
given set of parameters shows limit-cycle behaviour, we calculate 
the time-averaged value of L rad for an integer number of complete 
limit cycles. 

For a 10 M Q black hole at low accretion rates, we find as ex- 
pected < M„b s >= M. With increasing accretion rate, however, the 
advection comes into play and the "observationally calculated" ac- 
cretion rate deviates from the actual M. For a 10 M Q black hole 
this relation is practically independent of a. This deviation may be 
related to the deviations fro m the L rat j oc 7*^ relationship recently re- 
ported bv lKubota & Makishimjj2004l) and lAbe et ai]j2005l) . They 
find the so called 'apparently standard' regime in observations of 
XTE J1550-564 and 4U 1630-47, where the disc luminosity is pro- 
portional to 7? and attribute this to effects of the radial advection 
of energy. 

For the 10 6 M Q black hole (Fig.Q. There are departures from 
< M„bs >= M . While the disc is being fed with an accretion rate of 
M = 0.1M Edd , < M obs > is only 44 % of the actual accretion rate, 
while the time average for the low-M state only is as low as 1.5 % 
of the actual accretion rate. This is a result of the strong outburst 
behaviour of these discs compared to the low-mass case. Most of 
the mass previously stored in the outer disc is pushed through the 
inner disc in the high-M state. Then, however, most of the energy 
created by viscous dissipation is advected into the black hole. Since 
the outbursts are short compared to the complete limit cycle, the 
efficiency is fairly low. 
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4.4 Power density spectra 

We calculate lightcurves at equidistant time points. Suppose we 
have a lightcurve covering the time T with N points. We calculate 
power density spectra using the canonical normalisation to get the 
Power P v in units of (rms/mean) 2 /Hz. Integrating P v dv gives the in- 
tegrated fractional rms. We compute the FFT of the lightcurve and 
then define P v by <Lewin et aiTl9 88) 



2|o£ 
L 2 T 



(46) 



where L is the mean luminosity of the lightcurve and a v the non- 
normalised Fourier coefficient. The P v 's are binned into logarithmi- 
cally equidistant points. For the FFT we use the routines from the 
FFTW library 2 . 

The integrated fractional rms is given in percent 



100% 



PyClV 



(47) 



Using Parseval's theorem, it can be shown, that r indeed cor- 
responds to the rms/mean of the luminosity fluctuations 



100% 



\ 



(48) 



where L, is the luminosity given at equidistant time intervals and L 
the time-average of the total light curve. Unless otherwise stated, 
the last formula is used to calculate the integrated fractional rms. 

4.5 Period of limit cycle 

A sample lightcurve and power spectrum for a 10 M Q black hole is 
shown in Fig. [6] We binned the FFT data in logarithmically equis- 
paced bins of width Av/v = 0.005. The disc undergoes limit cycles 
on a timescale of roughly 600 s. Note that besides the fundamental 
frequency and their higher harmonics there are more peaks present. 
These arise from the asymmetry of the lightcurve. 

The luminosity of the disc is given in terms of M \, s /M, where 
M is the external/average accretion rate and M obs the "observed" 
accretion rate. This ratio is equivalent to the ratio between the ac- 
tual luminosity and the luminosity of a stationary disc where no 
advection is taken into account. 



5 COMPARISON WITH PREVIOUS WORK 

In this Section we briefly compare our models with previous and 
recent work in the field. We restrict ourselves to optically thick disc 
models. 

There is no absolutely correct way of doing ID accretion disc 
models in the one-zone approximation, and in this paper we have 
adopted the simplest approach, by ignoring all the details of the 
vertical structure. To do it properly clearly requires doing two- 
dimensional hydrodynamics. In Appendix IaI we derive the LHS of 
the energy equation <1 H explicitly from first principles. We work 
in terms of the total internal energy and the mass of a disc annulus. 
All other variables can be expressed in terms of mass and internal 
energy of a disc annulus. 

In early work on black hole accretion discs, iHonma et alJ 

ll99ll) presented disc models with an energy equation containing 

2 http://www.fftw.org/ 
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Figure 7. Lightcurve for an accretion disc around a 10 M Q black hole 
accreting at 0. 1 Me<m ■ & = 0.1. The metallicity is Z = 0.1, the number of 
grid points has been set to N = 250. Note that the lightcurve resembles the 
instabilities present in Fig. [3] 



the evolution of the total energy. They apply correction factors to 
try to take account of details of the vertical disc structure. They 
mainly consider M = 10M o and a = 0. 1 but make use of a modified 
"a-/'"-description by multiplying the n^-component of the viscous 
stress-tensor with a factor f? p . (Here f} P is the ratio of gas pressure 
related to total pressure, so that q = corresponds to what we use 
here, and q = 1 corresponds to the n/>-component being propor- 
tional to aP gas .) For q > 0.5 they find stable disc solutions while 
for q < 0.5 they find limit-cycle solutions of increasing amplitude. 
Their q = results are comp atible with ours. 

I n a series of papers ISzuszkie wicz & Millet] 1 1 9971 1 1998L 
1200 ll) take a similar approach to iHonrffitTHriTI 19911) except that 
they use different numerical factors to try to take account of the 
details of the vertical disc structure. For this reason, while their 
results are in line with those of Honm a et alJ ll99ll) and those pre- 
sented here, they differ in some minor details. In particular, they 
find unstable solutions for a 10 M Q black hole and a = 10~ 3 for 
0.09 < L EM < 1 which is in line with our models except the high-M 
end, where we find stability at 1 Lem. They subsequently develop 
a scheme where they take account of the possibility that the disc 
might become opticall y thin and there is non-Keplerian rotation. 

iNavakshin et alJ 12000) present a model with the specific ap- 
plication to GRS1915+105. They additionally include thermal and 
radiation diffusion in radial direction. For the simulations presented 
in this paper, however, these terms never become important. They 
might become comparable to the terms containing the radial advec- 
tion of the internal energy, but this only will result in factors of two 
at most. In the context of working in the one-zone-approximation, 
the omission of these terms in our simulations seems to be justi- 
fied. They investigate models with different radially varying values 
for the viscosity parameter a. They furthermore try to account for 
additional cooling in a corona and include some flickering. Their 
flickering is a random ad-hoc modulation of the efficiency of radi- 
ation coming from the inner parts of the disk. Lightcurves similar 
to GRS 1915+ 105 can be produced. No further diagnostics of the 
flick ering process (i.e. p ower density spectrum) are tested there. 

Janiu k*etall 120021) consider a similar approach to the one pre- 
sented in this paper. They write the change in specific heat in terms 
of temperature and density (i.e. surface density and scale height, 
setting the radial advection of the scale height to zero). They apply 
yet another set of correction factors to try to take account of the 
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Figure 6. (left) Lightcurve (detail) and (right) Power density spectrum for an accretion disc around a 10 M Q black hole accreting at 0.5 Mem- ce = 0.1. The 
FFT data are binned in logarithmically equispaced bins of width Av/v = 0.005. 



vertical structure . Their results are in line with those obtained by 
ISzuszkiewicz & Milled Jl997lll998ll200ll) . 

In conclusion, these various approaches differ mainly in the 
nature and actual numerical values of the correction factors (being 
in the range 0.2 - 16) applied when trying to take some account 
of the vertical disc structure. There seems to be no unanimity on 
the values that should be applied. In this paper we have taken the 
simple approach of setting all these factors to unity. 

Some of the work reviewed here iNava kshin et alJ feoOO: 
Janiu k et alj|2002h also include prescriptions for energy loss to a 
wind. While there is agreement that mass loss does not influence 
the results, energy loss (cooling) does influence the disk and sta- 
bilises it. Energy loss usually is only parametrised with a fraction 
of the radiative losses taken away to the corona or the wind. This 
fraction is either set constant or depends on the accretion rate (e.g. 
Ijaniuk et alj|2002l) . Without the energy loss to a wind, our results 
so far are in agreement with the results in literature. 

In the next section we shall introduce our model for the flick- 
ering which contains a self-consistent coupling between the mass, 
energy and angular momentum loss of the disc to a wind coupled 
to the physical model for the flickering. 



6 THE MODEL FOR THE FLICKERING 
6.1 The basic mechanism 

We now address the means of introducing a physical mechanism to 
produce the observed short-timescale variability in black hole ac- 
cretion discs. As we mentioned in the Introduction the main prob- 
lem is to find a mechanism which gives the right amplitude and the 

right times cale. 

iLvubarskill 1 1997ft showed that if the viscosity parameter a 
changes in a spatially uncorrelated manner on the local viscous 
timescale at each radius, then the characteristic flickering spectrum 
and the long timescales can be produced, but was unable to suggest 
a phys ical mecha nism which would p roduce this result. lKing et alJ 
(2004), following iLivio et alJ feoojli . put forward a solution. The 
main problem is that local variations in a are correlated with small 
variations in the disc's internal magnetic field which is generated 
by a local disc dynamo. This operates characteristically on about 
the local dynamical timescale. They envisage the disc dynamo op- 
erating as a set of essentially independent dynamo cells with char- 



acteristic radial width comparable to the local disc thickness. They 
then propose that the dynamo mechanism also gives rise to a small 
poloidal field in each cell. Most of the time this local poloidal flux 
is randomly aligned from cell to cell and so has no net effect. But 
from time to time this poloidal field magnetic field is sufficiently 
aligned from cell to cell that it gives rise to a large scale magnetic 
field which is able to generate an outflow (wind or jet) which re- 
moves angular momentum from the disc. Since the local field in 
each dynamo cell changes on about the local dynamical timescale, 
this overall alignment only happens on a large multiple of the dy- 
namical timescale and thus can be comparable or even longer than 
the viscous timescale. 

We now apply these ideas to the thermal disc structure de- 
veloped above. We stress here the additional concepts required in 
applying this model to a realistic disc. More details of the underly- 
ing ideas can be found in lKing et alJ 120041) . In the following sub- 
sections we describe the details of this model and the equations 
involved. 



6.2 Evolution of the poloidal magnetic field 

We consider a poloidal magnetic field B- which evolves according 
to the induction equation 

dB- 15/ ,dB,\ 1 d 

We allow the field to diffuse radially with an effective magnetic 
diffusivity if and allow for radial advection at velocity U R which 
is due to a magnetic wind torque. S B stands for a source term for 
the poloidal field which is assumed to be generated by the local 
magnetic dynamo. At the outer disc boundary we take B- = 0, 
allowing magnetic flux to diffuse outwards, but preventing inward 
flux advection. 



6.3 The source of poloidal field, S B 

The poloidal field B- is assumed to be generated by the disc dy- 
namo which generates the main disc viscosity. The magnetic field 
in the disc, partitioned into disc annuli of radial extent AR ss H, is 
assumed to change independently in each dynamo cell in a stochas- 
tic fashion on the characteristic dynamo time scale 



(50) 



10 M. Mayer & J. E. Pr ingle 




-0.5 0.5 
log (v x d ) 

Figure 8. The power spectrum is shown for a sample time-series of a 
Markoff chain series (see eq. 1531 for &[ = 0.5. This series is used to mimic 
the effects of locally acting magnetic dynamoes. In this figure we show 
the power spectrum for a dynamo acting on a timescale tj. To account for 
timescales smaller than tj, we linearly interpolated the Markoff chain. The 
frequency is normalised to tj. The peak at the dimensionless frequency 
log(vr</) = is clearly visible including higher harmonics. Most of the 
power is released on slightly longer timescales (approximately 2.5tj) com- 
ing from the broad peak at log(vTrf) ss -0.4. 



which we take to b e a mult iple of the loc al dyna mical time scale 
(we take k A = 10, cf. lTout & Pringlell99llstone et alJl99r3) 

Each annulus of the disc is assumed to generate a suf- 
ficiently large internal tangled magnetic field of strength B disc 
which ge nerates the effec tive disc viscosity with a parameter 
JShakura&Sunvaevll973lt 



g disc 

AnP' 



(51) 



Although we assume that each dynamo cell generates a local 
poloidal field B-, we expect the magnitude of the poloidal com- 
ponent to be small. We define a parameter B s so that the maximum 
value of the poloidal field, B z mas is given by 

Bzjnax /^s^discs (52) 

where 6 disc is given by Equation l51l We expect B s <k 1 . 

We then model the stochastic nature of the local poloidal field 
in the following way. At each radius we define a set of times, t„ = 
riTd, n = 1, 2, ... as integral multiples of the local dynamo timescale 
at which the field changes. At each of these times we generate a 
random number, «(?„), according the the Markoff process 



«fc,+i) = -£*l"fe) + efo) • 



(53) 



where e(f„) is a random variable of zero mean and unit variance, 
and we take a { = 0.5. The new value of the magnetic field strength 
in the disc annulus then is given by 



B z (t„+i) - u(t n+ \)B-„ 



(54) 



We show a power spectrum of a sample time-series for u(t„) 
for a magnetic dynamo acting locally on the dynamo time scale r d 
in Fig.|8| For the purposes of the Figure, in order to resolve shorter 
timescales than the dynamo timescale, we applied a linear inter- 
polation between the timesteps. We note that the power spectrum 
peaks at a period of a few times tj. 

As mentioned above, in order to ensure the consistency of the 
model, we need to make sure that the poloidal field B- is always 



smaller than the disc magnetic field. This implies Bj, <k 1, so that 
energetically the poloidal field is negligible. Thus this assumption 
furthermore allows u s to neglect the energy generation/loss due to 
our dynamo process. Tout & Pringle 1 1996) estimate this fraction 
to be ~ H/R. 



6.4 Magnetic torque and the wind 

We assume that angular momentum loss occurs due to a magnetic 
wind/jet when the poloidal field is of sufficiently large scale. This is 
likely to occur when the poloidal fields generated by the individual 
dynamo cells are, by chance, spatially correlated over a radial ex- 
tent of the order of R. To measure the degree of spatial correlation 
in the simulation, we define the radial averages 



<S : > 



B: 



(55) 



(56) 



La- RdR 

We choose the total interval over which we measure the cor- 
relation as equal to R, that is we take 

A" + A + = R. (57) 

We then choose R to be the geometric mean of R - A_ and R + A + . 
Thus we choose the dimensionless constant a so that R- A" = R/a, 
R + A + =aR with a > 1. This then implies that a = (1 + V5)/2, the 
golden ratio. 

Thus we obtain 



(B : ) 



f R+A+ B-RdR 

JR-A- 

R 2 (a 2 -l/a 2 )' 
f R+ f B 2 RdR 
R 2 (a 2 -l/a 2 ) 



where we note that a 2 - I /a 2 = VB". 
We now define the quantity 



Q 



(B z f 



(B 2 



(58) 



(59) 



(60) 



which represents the degree of spatial correlation of the poloidal 
magnetic field in neighbouring disc annuli. We note that < Q < 1, 
where 2=1 represents a fully coherent magnetic field, while for 
2 = the magnetic field is fully randomly ordered. 

The poloidal magnetic fie ld B 7 produces a torqu e 7" mag on a 
disc annulus, width AR, of size fcjvio & Pringle 1992) 



-AnR 1 



An 



AR. 



(61) 



Considering the change of angular momentum in the disc an- 
nulus, this can be combined to give 

2nRARj t (l/? 2 n K ) + A (2tiRI,v r R 2 £1 k ) = AG + r mag , (62) 

where the terms on the LHS describe the local change of angular 
momentum and the radial advection of angular momentum from 
neighbouring annuli, while the terms on the RHS describe the vis- 
cous and magnetic torques. Note that while angular momentum 
transport due to radial advection and the viscous torque for the 
whole disc reduce to the values at the inner and outer boundary, 
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the magnetic torque produces a local angular momentum loss term 
across the disc. 

Given this additional term (r mag ) in the angular momentum, 
equation J3) now becomes 

d i , \ 1 d i j \ 1 8G l(B z ) 2 

— (£# 2 n K ) + [RZ.v R R 2 n K ) = 2R 1 

<9r v > RdR y > 2nRdR 



4tt 



(63) 



In contrast to lKing et alJfcOO^h . we now take explicit account 
of the fact that the wind/jet, which removes angular momentum 
from the disc, may also remove a significant amount of mass. We 
assume that the magnetic outflow at radius R co-rotates with the 
disc along each field line, out to some Alfven radius _Ra(> R)- Then 
the local mass loss rate, Lj, is related to the torque by 



IkRARLzRICIk = T m 



which implies that 



Lv = - 



2(R/R A f j(B z ) 



CIkR 



4k 



(64) 



(65) 



The loss term leads to a modification of the continuity equa- 
tion in the form 



<3£ 1 d , ^ s , 
Tr + -RTR (KLVR) - U - 



(66) 




Eqns. l66l and 1631 now can be combined using the functional 
form of the viscous torque (eq.[4j to give 

— = Vfl — ( v t E Vr) 

dt RdR\ dR^ > 



(67) 



For R = R A the enhancement of angular momentum due to the ex- 
tra torque r m;lg is balanced by the loss of angular momentum taken 
away by the wind since it is taken away with the same lever arm 
(see eq. 1641 and the poloidal magnetic field influences the evolu- 
tion only through the mass loss. Note that iKing et alj 120041) es- 
sentially assumed that = °°, and thus that mass loss is not im- 
portant in the model presented there. However we retain the ratio 
Ra/R as a parameter and allow for mass loss, since analytical esti- 
mates IPudritz & Normarjfj 986 ) give Ra/R = 10 while numerical 
simulations lOuved & Pudrita l 19971) suggest R A /R ~ 1.5. Unless 
otherwise stated we choose Ra/R = 3 as the fiducial value. 

The mass lost also removes the internal energy of the gas at a 

rate 



(68) 



and thus we also require a revised energy equation (see eq. lilt 
which is in the form 

dq de d(2nRHv R ) 

— = e + v R — + 2nRARPH + P— — 

* dR dR (69) 

= 2nRAR(Q + + Q w -Q). 

We now have three equations (eq ns . l67l 149 1 and |69 \ describ- 
ing the evolution of the disc. These are integrated in time using 
a one-step Euler scheme. Advection is treated in a first-order, up- 
wind donor cell procedure. We note that the main observational 
output from the disc is the radiation L mi coming from the disc (see 
eq.|42j. 



6.5 Magnetic diffusivity 

Following lKing et alJ 120041) . we take the effective magnetic diffu- 
sivity if in eq.|39]to be 



77* = v t max(l, QR/H), 



(70) 



where we assume that the Prandtl number is of th e order unity fa = 
v t ), but allow it to be enh anced by a factor R/H Ivan Ballegooiien 
ll989tlLubowetalJll994l) and include large scale effects by mul- 
tiplyi ng it by the coherency factor Q. For further discussion see 
lKingetalJh004 . 



6.6 The interplay between the timescales for the flickering 

In this section we quantify the timescale argument for the appear- 
ance of flickering. As mentioned previously, fluctuations propagate 
through the disc, if the viscous timescale is shorter than the mag- 
netic timescale, i.e, the timescale where we expect sufficient align- 
ment of the poloidal magnetic field in neighbouring cells. 
The viscous timescale is given bv jPringlell98l[) 

1 



."Tdyn , (71) 

a (H/R) 

where Td yn = 1/Q.k stands for the dynamical timescale. H/R is 
the opening angle of the disc. The magnetic timescale, i.e. the 
timescale, on which we expect neighbouring field lines to be suffi- 
ciently aligned, is given bv jLivio et all200i 

Thus the condition for propagation is 



(72) 



ag/ T visc ^ 1 



(73) 



which translates into 



2 m k d a (H/R) 2 > 1 



(74) 

We plot the ratio r T in Fig.|5|for different values of a. The value 
of k d is fixed to 10 (see Sect. 16.31 . The ratio is smallest for 
(H/R) mm « 0.35. For a = 0. 1 fluctuations formally propagate only 
for H/R > 0.5 or H/R < 0.25 only, while for a = 1 fluctuations can 
propagate for all H/R. In contrast to that, for a = 0.01 fluctuations 
only propagate for H/R < 0.1. 

This is in line with findings of IChurazov et alJ 1200 lh who 
state that fluctuations close to the dynamical timescale are effec- 
tively damped in geometrically thin discs but can survive in geo- 
metrically thick discs since then the fluctuation timescale becomes 
comparable to the viscous timescale. In our model fluctuations for- 
mally survive for lower H/R ratios as well. Then, however, the 
magnetic timescale is so long that alignment becomes extremely 
rare. From the above analysis this probably puts a lower limit on 
the viscosity parameter or of a ^ 0.1 since most of the variability 
is expected to come from a optically thin, but geometrically thick 
disc. 



7 OBSERVATIONS OF THE FLICKERING IN THE 

THERMAL COMPONENT OF THE X-RAY SPECTRUM 

Before describing the numerical simulations we first consider the 
likely observational scenario to which they are relevant. In gener- 
alising the work of IKing et alJ 120041) to more physically realistic 
discs, we have initially confined ourselves to considering the stan- 
dard accretion disc configuration, with the addition of heat advec- 
tion. In general this is assumed to correspond in X-ray binaries to 
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Figure 9. Ratio r T of the magnetic and viscous timescale in dependence 
of H/R (see eq. 1741 for different values of the viscosity parameter a and 
kj = 10. Only for r T > 1 fluctuations can propagate. 

the thermally dominant (TD) or high/soft state. However, from an 
observational point of view, the variability in the TD or high/soft 
state is much less than in the LH state. 

Usually the main components in the spectrum of X-Ray bina- 
ries caii_be_well_^t_by_j_black-body and a power-law component 
(cf. iMcClintock & Remillard J2004h ). The black-body or thermal 
component for galactic X-Ray binaries appears at energies around 
1 keV, while the power-law component extends to higher energies. 
The components are thought to represent the radiation from the op- 
tically thick disc and a optically thin corona, respectively. Thus, 
since in the current paper, we only attempt to describe the time- 
dependent behaviour of the thermal disc component, i.e. the opti- 
cally thick disc, we need to restrict ourselves to observations which 
explicitly concentrate on measuring the flickering in the thermal 
disc component of the X-R ay spectrum. 

Chura zov et alj |2001 ) show a power spectrum density (PSD) 
of Cyg X-l in the high-soft state for energies of 6-13 keV, and state, 
that the amplitude in the softer bands is much smaller due to the 
influence of the black-body emission from the disc. They argue that 
their data is consistent with essentially all of the variability being 
associated with a power-law component and not with black-body 
emissi on (i.e. disc compone nt), but do not quote a formal limit. 

Mivamoto et al. 1 1994) examine Ginga observations of Nova 
Mus 1991 (GS 1124-683). They disentangle the variability in the 
disc and power-law component by calculating the PSD for differ- 
ent observations which have different fractional contributions from 
the power-law component to the total flux. For a power-law com- 
ponent fraction above approximately 10 per cent, the PSD at 0.3 
Hz is an increasing function of power law fraction. While for lower 
values it is roughly constant (1.5-10" 5 Hz" 1 ). At these low values, 
the slope of the PDS is about -0.7 but due to the error bars it still 
could be fit with -0.5. Nowak et al. 1 1999) show that for counting 
noise dominated time-series the expected PDS slope is -0.5. Thus 
we conclude the data are consistent with a non-detection of flicker- 
ing in the disc component. By using the value of the PDS at 0.3 Hz 
only we estimate an upper limit to the rms variability of r < 0.2 per 
cent. 

iNowak et alJfeOOll) consider LMC X-l and LMC X-3. These 
are persistent black hole sources (i.e. not transients), and generally 
show spectra dominated by a soft, thermal component. They detect 



no variability in LMC X-3 down to a level of about 0.8 per cent. 
They do detect variability in LMC X-l of 7 per cent throughout 
the spectral range considered (0-9 keV), but since for LMC X-l the 
power law component is much stronger that in LMC X-3, this could 
just be variability of the power law component. 

Observations of the transient 4U 1543-47 during its 2002 out- 
burst are presented in IPark et all2004h . During most of the decay 
the spectrum is soft and dominated by black body emission. The 
PDS shape in this phase is roughly 1/v and the variability is about 
r = 1 per cent. However the fraction of black body emission never 
exceeds 80 per cent, i.e. the power law always co ntributes more 
than 2 per cent. Comparison with the results of iMivamoto et all 
1 1994) shows, that for a black-body flux fraction of 20 per cent, if 
the relation shown there is representative, the fractional rms must 
be r = 1 per cent. 

iKalemci et alJ<2004l) report on the light curves of a number of 
transients observed with RXTE. They again find strong variability 
in the low/hard state, while in the high/soft state they only can give 
upper limits, since they are counting noise dominated. Depending 
on the source, and the quality of the data, these limits are in the 
range 1 to 8 percent. 

To quantify the relation between the flux fraction a nd the 
rms variability, we hav e taken data from K alemci et alj 120041) and 
IMivamoto et alj 1 19941) and plot the rms versus the black body flux 
fraction in Fig. llQI The fluxes are calculated in the 3-25 keV band. 
Since lMivamoto et alj 1 19941) only give the normalised power den- 
sity P{v) at frequency v = 0.3 Hz, we estimate the integrated rms r 
by using r = 2 yfvP~ v . The factor 2 is introduced to account for the 
dv/v in the definition of the integrated rms (cfeqj47l. The black- 
body fraction for the LMC X- 1 and X-3 data of lNowak et alj <200ll) 
has been calculated using their black-body + powerlaw fits. While 
there is some scatter in the data, the global trend is that for black- 
body flux fractions lower than 5 per cent the rms variability is about 
25 per cent, while it drops beyond detectable limits for larger val- 
ues of the black body fraction. The transition is not very sharp, but 
occurs around a black body fraction in the range 3 to 20 per cent. 
The transition might be sharper than this but we are only using data 
in the 3-25 keV range where the black-body contribution is small 
compared to the power-law compone nt, and varies fro m source to 
source. In addition, we note that Zdziarski et al. 1 2005) carry out an 
energy-dependent analysis of the fractional rms for GRS 1915+105 
and come to similar conclusions. 

To summarise, it is clear from the observations that the rms 
variability of the black-body component in the X-Ray spectra of 
black hole sources is small. Indeed there are no clear detections, 
but only upper limits. As a representative value we take r = 1 per 
cent as an upper limit for the variability of the thermal component. 



8 RESULTS 

To demonstrate the effect of the stochastic magnetic wind/jet on our 
disc models model with the physical input so far, we ran some mod- 
els to assess the influence of the parameters on the variability. For 
all models, we fix the mass of the black hole to be 10 M Q , the vis- 
cosity parameter a = 0. 1 (a ppropriate for standard optically thick 
discs, see lLewin et alJl 997) and take for the dynamo timescale (cf. 
eq.|50j k d = 10. 

First we consider a black hole accreting at half of the Edding- 
ton rate, M = 0.5M E(M . We take the Alfven radius at each radius 
to be a constant value of R A /R = 3. From Figure 4, in the absence 
of magnetic flickering (i.e. with /?s = 0) the disc is unstable and 
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Figure 11. Sample lightcurves for a 10 M G mass black hole accreting at 0.5 Mryd f° r different values of /?s, the strength of the poloidal field compared to 
the disc magnetic field. The ratio of the Alfven radius Ra to the radial distance R is R&/R = 3, the viscosity parameter a = 0.1, The intervals indicate the 
time-segments shown in Fig. ll2l used to calculate the power spectra shown in Fig. l 1 31 
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Figure 10. Fractional rms variability r for different black-body flux frac- 
tions i n the 3-25 keV band for different objects. All data from Kalem ci et alJ 
1 2004) except the dat a for GS 1124-683, which have been taken fro m 
Miya moto et"a?]<1994h and LMC X-1 and X-3 from lNowak et alJ l200ll) . 



shows limit-cycle behaviour. Fig. ll ll shows sample lightcurves for 
different values of /?s. The limit cycle is still present, i.e. the disc 
oscillates between a high M a \, s and low M a \, s state. It is best evi- 
dent from the /3 S = 1 case that while in the high M obs state there 
is some variability, in the low M b s state there is hardly any sign 
of variability. This is a result of the smaller H/R ratio in the in- 
ner disc during the low M obs state. Then the alignment timescale 
(proportional to the factor 2 R/H times the dynamical timescale, see 
iLivio et alJl2003l) is much longer than the viscous timescale (pro- 
portional to (R/H) 2 times the dynamical timescale): Fluctuations 
formally could survive, but since alignment occurs on extremely 
long timescales only, they are un likely to produce flickering (see 
Section and iKing et aljEool) . In the high-M state, however, 
both timescales are at least comparable and so fluctuations propa- 
gate to the inner disc without being smeared out by the viscosity. 
It is evident from Fig. 1111 that a non-zero /?s influences the 



strength of the variability. In Fig. ^| we show the influence of the 
parameter /?s on the power density spectra (PDS), while the inte- 
grated fractional rms r (see eq. !47> is shown in Fig. 1141 

The time segments used for the calculation of the PDS are 
indicated in Fig. II II Throughout this section, we only take seg- 
ments from the lightcurve, where the disc was in the high-M state. 
Since for red-noise the lightcurve is only a stochastic realisation 
of the underlying process, the calculated PDS fluctuates randomly 
around the "true" spectrum (e.g. lUttlev et alj|2002l) . To overcome 
this difficulty, we divided each segment in 10 equally spaced sub- 
segments, calculated the PDS of each subsegments and took the av- 
erage of these. We further normalised each subsegment to remove 
linear trends in the lightcurve. 

We see from Fig. 1 141 that the rms variability is r < 1 %, for 
/3s < 0.2 and that r strongly increases for larger /3s and reaching 10- 
20 % for /?s = 1. The near quadratic increase of r with /?s reflects 
the fact, that the poloidal magnetic field influences the disc only 
through the torque which is proportional to which in turn is 
proportional to p\P, i.e. a fraction f}\ of the viscous torque (see 
eans. l6ni52ll5Tl andl22l 

Given the observational constraints, discussed in Section 7, it 
is evident that we need to take (3$ to be smaller than around 0. 1 - 
0.2. Thus in this model we require the energy density in the poloidal 
field component to be at most only a few percent of the energy 
density in the magnetic field generated by MRI in the disc. 

Next we consider the influence of the accretion rate on the 
variability for constant y6s- Fig. ll5l shows the PDS for two values of 
/3s for different accretion rates. While the integrated rms variability 
r does not depend on the accretion rate but on a, the shape of the 
PDS for different accretion rates is slightly different. 

We also briefly investigate the behaviour which occurs if /3s, 
the ratio of poloidal and tangled magnetic field strengths, is as- 
sumed to_var^_wjth///i? ; _For example, in their simple dynamo 
model, Tout & Pringle ( 1 996 ) suggest that /?s might be proportional 
to H/R. If so, since H/R in the simulations described is at most 
H/R ~ 1/5, we can put a constraint on the constant of proportion- 
ality to be smaller than unity. 

In Figures [T6landfT7l we show the results of simulations with 
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Figure 12. Lightcurve segments for the cases shown in Fig. 1111 for dif- 
ferent /?s ■ These are used to calculate the PDS in Fig. 1131 M„|, s is given 
in arbitrary units, but the amplitude still has got some meaning. Note that 
the overall shape of the lightcurve does not change considerably, while the 
range in M b s changes significantly. The longterm trend in the data has been 
removed. 



the same parameters as before but with p s = 0.5H/R and p s = H/R 
for varying accretion rate. For the 1 Mem cases there is a significant 
decline of the PDS for small frequencies. This decline is a result of 
a declining H/R in the outer disc which shortens the amplitude of 
flares released from further out (i.e. longer timescales) through the 
///^-dependence of /?s- For the lower accretion rates this decline is 
far less obvious, but then the time the disc is spending in the high- 
M state is becoming shorter and so does the time available for the 
FFT and averaging. The integrated fractional rms in for the same 
accretion rate changes by a factor of about 4 when comparing the 
cases of y6 s = 0.5H/R and f3 s = H/R consistent with the result of 
Fie. ll4l where we show the quadratic dependence of r with /3s • 

Finally we explored the influence of different values of Ra/R 
on the variability. In Figure IT81 we show the influence for a 10 
M black hole accreting at 0.5 M Edd with /3 S = 0.25. The sys- 
tem changes with Ra/R increasing from a mass loss to a mag- 
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Figure 13. Power spectra for a 10 M G model solar mass black hole, accret- 
ing at 0.5 M^M f° r different values of /?s for the lightcurve segments shown 
inFig.IHlff = 0.1 
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Figure 14. Integrated fractional rms r for a 10 M model solar mass black 
hole, accreting at 0.5 Meaa in dependence of /3s, the strength of the poloidal 
field compared to the disc magnetic field (cf. Fig. 1131 . The viscosity param- 
eter is a = 0. 1 . 

netic torque dominated regime. For the parameters chosen, the in- 
tegrated rms-variability r is increasing with Ra/R and saturates 
for Ra/R > 3. The power spectral shape then is indistinguishable 
within the error bars. Thus the assumption of iKing et alJ 1 2004) 
who set Ra/R = °° is justified as long as the "real" Ra/R > 3. 



9 CONCLUSIONS 

We have generalised the model o f lKing etalJ l2004) for variability 
in black hole accretion discs by including proper consideration of 
the local disc structure. We take the disc properties to depend only 
on radius R, and so, in common with other authors, make a local 
one-zone approximation for the disc structure. 

In the absence of the addition of a stochastic magnetic dy- 
namo, which drives the flickering, we mainly reproduce results 
similar theoretical work has shown. We agree that taking account 
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Figure 18. (left) Power spectra for different values of R\/R, the ratio of the Alfven radius R\ to the radial distance R for a 10 M Q black hole accreting at 
0.5 Mem with /3s = 0.25, the ratio of the poloidal field compared to the disc magnetic field, (right) Integrated rms-variability r as a function of R\/R. For all 
Rx/R > 3 the power spectra and the value of r is essentially indistinguishable. 
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Figure 15. Power spectra for a 10 M Q model solar mass black hole, accret- 
ing at different accretion rates for different /?s = 0.125 and 0.5. a = 0.1. 




Figure 16. Power spectra for a 10 M Q model solar mass black hole, accret- 
ing at different rates for/3 s = 0.5 (H/R). a = 0.1. 



to advective heat flow is of crucial importance in determining the 
structure of the inner regions, but find that the assumption of strictly 
Keplerian angular velocities, and the neglect of radial pressure gra- 
dients, are generally reasonable assumptions. 

When a stochastic dynamo is included we find behaviour sim- 
ilar to that described by K ing et all i2004V However, because the 
degree of variability of the thermal disc is observed to be small (less 
than around 1 per cent), if it is detected at all, we can only draw 
limited conclusions. Our main finding is that for consistency we re- 
quire that the strength of the poloidal field generated by any radial 
disc dynamo cell must be at least an order of magnitude smaller 
than the field generated by the dynamo within the disc. The simula- 
tions reported here however are of help to understand the influence 
of the parameters on the variability. 

In this paper we have restricted ourselves to modelling just the 
thermal disc component, and so have had to restrict our attention to 
the high/soft or thermally dominant state of black hole binaries. In 
future work we shall begin to include consideration of a hot corona, 
in addition to the cooler thermal disc, so that we can begin to model 



the more highly variable low/so ft state of th e black hole binaries, 
as well as the X-rays from AGN lUttl ev et all2005l) . 
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APPENDIX A: THE ENERGY EQUATION 
Al Derivation of the LHS 

The first law of thermodynamics is 

dQ = dE + Pd^j, (Al) 
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where dQ stands for the change in heat Q per unit mass, dE the 
change in internal energy E per unit mass, and Pd(l/p) the volume 
work. P is the pressure and the specific volume is 1 /p. 
The continuity equation in general form reads 



d(pV) = pdV + Vdp = 0, 
which is equivalent to 

dV _ dp 
~~~~p~- 

In total differentials using the continuity equation, the energy 
equation can be shown to be equivalent to 



(A2) 



(A3) 



d(pVQ) = d(pVE) + PdV. 



(A4) 



This relation implies that the change of total heat in a volume V is 
given by a change in the total internal energy and the expansion of 
the volume. 

In cylindrical coordinates, vertically averaged, we have 



Thus 



V = 2nRHAR (A5) 



d (tiRZARQ) = d (nRLARE) + 2nRARPdH. (A6) 



Replacing the specific heat Q by the total heat content q in a 
semidisc annulus, 



q = nRZARQ, 



(A7) 



and the specific internal energy E by the total internal energy e in a 
semidisc annulus, 

e = nRZARE , (A8) 

we can write 

dq = de + Pd (2ttRHAR) . (A9) 
Time-dependently, the LHS of the energy equation is written 

(A10) 



dq _ de p d(2nRHAR) 
dt dt dt 



Taking into account advection, with radial velocity u R , 

dq de . d(2nRHu R ) 

— = e + u R — + 2nRARPH + P— — . (All) 

dt 8R 8R 

In discredited form, 

— = e + A {nRusZE) + 2nRARPH + PA (2nRu R H ai ) , (A12) 
dt 

where // ad is the adverted value of H, and this term corresponds to 
the adverted volume contributing to the 'PdV work. 
Now, given the specific internal energy E, where 

3 kT Act , 

E=~ + —T\ (A13) 

I fiiiip cp 

and given hydrostatic equilibrium in the form 



kT Act. 



M 



P + -^~ T * = G- — r; 

fim p 3c ApR i 



(A14) 



p and T can now be be computed, for given values of E and S. The 
variation of H in terms of E and £ is given by using the definition 
of the surface density (Equation^ 

f = f-^. (A15) 
H L p 

Using both the hydrostatic equilibrium, and the definition of the 



internal energy E, we can express the variation in p in terms of 
variations in the surface density S and specific internal energy E as 



dp 



-6i] dS, 4-3/3/. dE 



(A16) 



p 8 + fit - 7/7 S 8 + p P - 7/7 E 

With the combination of the last two equations we can eliminate 
dp /p to obtain 



dH Pp-t] dL A-3p P dE 

IT ~ 8 +p P - 7/7 T + 8 +p P - 7/7 ~~E' 
APp-A-qdL A- 3p P de 
S+Pp-lqY %+Pp-l^~e , 
API -npp + S dE -3P 2 p + Wp -& de 
Pi + ttpp - 16 S + pi + Ep~ P - 16 T' 



where we use 



and 



(A17) 
(A 18) 

(A19) 
(A20) 

(A21) 
(A22) 



as the fractional contribution of the gas component to the specific 
internal energy and pressure, respectively. These two variables are 
related by 

The specific internal energy and pressure are related by 

E =-p( 3 - 3 2^)- (A24) 

Equation lA19l multiplied by H and taking the time-derivative 
leads to (using Equations <A23> and Q) 



Ap 2 p -Upp + 8 i -3p 2 p +l0p P -i 
H = + 



Pi + l3p P - 16 2p pi + Y3p P - 16 2nRpARE 



(A25) 



C E C e 

The coefficients C e and — Cj are equal to I for pp —* and 
Pp — » 1, while going through a minimum for p P ss 0.75 at values 
0.38 and 0.35, respectively. Then in both the limiting cases of p P — » 
and Pp — » 1, we find H oc ^feJT. ss •^Lcj JT, = c s , recovering the 
hydrostatic equilibrium. 

Using Equation <A25t in the energy equation <A12t . we can 
write 

— = e + A (7cRu R 2,E) + 2nRARPH + PA (2nRu R H^) 
dt 

I P \ P 

= e\l+C e — + A {nRitplLE) + nRAR—C^L + 2tcRPA (u R H. d6 ) 
\ P E I P 

= ell +C e — | + A\-ME.J + jiRAR-C^t. + -PA\M— |, 
\ pE) \2 L J p 2 \ p ad / 

where we applied M = 2ixLRu R in the last step. 



APPENDIX B: LOCAL STABILITY ANALYSIS 

The thermal stability analysis <Pringleil976l) . now including advec- 
tion shows thermal instability for 

faiog(2 + \ (dlog(Q- + Q aA ) 



dlogT 



dlogT 



> 0. 



(Bl) 



1 8 M. Mayer & J. E. Pr ingle 



Thermal instability occurs, if the increase(decrease) of temper- 
ature compared to the equilibrium state leads to a stronger in- 
crease(decrease) of the heating compared to the cooling while 
keeping hydrostatic equilibrium. 

Using the equilibrium condition (Q + = Q~ + gad) and 



we get the condition 



f a log Q + 
\ aiogr 



d log gac 

aiogr 



- (1 - e) 



d log Q- 
dlogT 



(B2) 



> 0. (B3) 



Along similar l i nes, viscous instability (e.g. 
Light man & Eardle^ 1 19741 : IPrineld Il98ll) only occurs as long 



<91ogM 
dlogE 



< 0. 



(B4) 



The disc is viscously unstable, if an increase(decrease) of accretion 
rate does lead to a lower/higher surface density. 

For the calculation of the criterion we use the hydrostatic equi- 
librium 



io g p-iog^5: 2 = / 1 (p,r,E) = o, 



(B5) 



and the stationary energy equation 

log Q + - log (Q- + g ad ) = f 2 (p, r,D = . (B6) 

With the contribution of gas and radiation to the total pressure 
(see ea. H6> we can write the total variation of the hydrostatic equi- 
librium and energy equation in terms of p, T and E, generalising the 
method of Maver & Duschlj 120051) for a non selfgravitating, opti- 
cally thick accretion disc. We neglect changes in Xa6 since Xad is 
expected to vary only very little. 

We take the variation of the hydrostatic equilibrium (f\ ) 







dp dT dL 
A—+B — + C — 

p T S 

and the stationary energy equation (/ 2 ) 

dp dT dL 
D— +E — + F — = 
p T S 

with respect to p, T and S. The coefficients A ... E are 



(B7) 



(B8) 



A = 

B = 

C = 

D = 

E = 

F = 
with 



aiogp 

a/i 
aiogr 

d logS 

aiogp 

dfi 
dlogT 

dfi 

a logs 



1+ySp (B9) 

4-3/3/. (BIO) 

-2 (Bll) 

A R (l-6 ad ) + ^ P -l)(l-2e ad ) (B12) 

(B R - 4) (1 - e ad ) + (4 - 3B P ) (1 - 2e ad ) (B13) 

2(1 -Qd) (B14) 



(Ar,Br): 



d\ogK R 
dlog(p,TY 

Thermal instability exists, if 

AE - BD > 



(B15) 



(B16) 



where we used the variation of the energy equation and expressed 
the density variation dp/p in terms of temperature dT/T using the 
variation of the hydrostatic equilibrium while keeping the surface 
density X=const. 

We have viscous instability for 

CE -BF AF - CD 

-(Bp -I) (4-3pV) + 1<0 (B17) 

^ AE-BD v AE - BD v ; 

where we both expressed the density and temperature variations in 
terms of the surface density variations and use this expression in 
the angular momentum equation 



M - 37rv t E = 



(B18) 



to get the instability criterion. 

In the radiation pressure (B P = 0) and Thomson scattering 
(Ar = Br = 0) dominated regime we get thermal instability only if 



4 - 12e ad > 



(B19) 



< 



(B20) 



and viscous instability as long as 

7fe ld - 1) 
1 - 3e ad 

Since the nominator of the LHS of the last expression is always 
negative (0 < e ad < 1), the condition of thermal instability and 
viscous instability are the same. 

To conclude, thermal and viscous instability only occurs, if 
fad < 5 • A radiation pressure dominated disc can be thermally and 
viscously unstable as long as advection only contributes less than 
one third of the local energy loss. This clearly shows the stabilising 
effect of energy advection on radiation pressure dominated accre- 
tion discs. 

For bound-free and free-free absorption (A R = 1, B R = — |), 
it can be shown that in this case the disc is stable with respect to 
thermal and viscous instability for all B P . 

It needs to be stated that these results are only based on a lo- 
cal stability analysis and it still needs to be shown by either time- 
dependent simulations or a global stability analysis when these un- 
stable modes are operational. 



